//package Powerflow;
//import BaseCalculate.*;

/**
 * 牛顿——拉夫逊法潮流计算程序
 * @author Administrator
 */
public class Powerflow {
	private static Plural[][] Y;//节点导纳矩阵
	private static Plural[] U  ;//各节点电压
	private static Plural[] v;//未设置平衡节点时，输出电压空值
	private static Mpqv[] PQU;//各节点的状态，有功，无功或电压值
	private static int n;//节点个数
	private static int s=0;//平衡节点号
	private static double[][] J;//雅克比矩阵
	private static double[] dPQU;//不平衡量
	private static double[] dU;//电压修正量
	/**
	 * 牛顿——拉夫逊法潮流计算主程序
	 * @param args 参数YY是一个Plural（复数）型的二维数组，存放节点导纳矩阵，
	               要求Y中的元素必须指向一个plural（复数）型的对象；
	 * @param args 参数UU是一个Plural（复数）型的一维数组，存放各节点电压
	 * @param args 参数PQU是一个Mpqv型的一维数组，
	               其中Mpqv型中包含各节点的状态，有功，无功或电压值
	 * @return 返回值是一个Plural（复数）型一维数组（各节点电压）
	 */
	public static Plural[] flow(Plural[][] YY,Plural[] UU,Mpqv[] PQU1) {
		Y=YY;
		PQU=PQU1;
		
		n=Y.length;
		U=new Plural[n];
		for(int i=0;i<n;i++){     
			U[i]=UU[i].clone();   //电压赋值
			/*if(PQU[i].mark==0){   //寻找平衡节点号
				this.s=i;
			}*/
		}
		/*if(s==0){                //未设置平衡节点时，输出电压空值
			return v;
		}*/
		int isCal=0;
		for(int i=0;i<n;i++){     
			if(PQU[i].mark==1){//寻找平衡节点号
				s=i;
				isCal++;
			}
		}
		if(isCal!=1){                //未设置平衡节点时，输出电压空值
			return null;
		}
	
		for(int i=0;i<10;i++){ 
				J=J();//求雅克比矩阵
				dPQU=dPQU();//
				dU=Mat.multi(Mat.inv(J),dPQU);//最多循环次数
				int sl=convergence(0.0001); //判断电压是否收敛
				U=New_U();
				if(sl==1){            //若收敛，则强制跳出循环
						break;
				}
		}
		return U;
	}
	/**
	 * 求节点注入电流的实部
	 * @param 参数i是一个int型的数（节点号）
	 * @return 返回值是一个double型的数（节点注入电流的实部）
	 */
	public static  double aii(int i){
			double aii=0;
			for(int k=0;k<n;k++){
					aii=aii+Y[i][k].a*U[k].a-Y[i][k].b*U[k].b;
			}
			return aii;
	}
	/**
	 * 节点注入电流的虚部
	 * @param 参数i是一个int型的数（节点号）
	 * @return 返回值是一个double型的数（节点注入电流的虚部）
	 */
	public static  double bii(int i){
			double bii=0;
			for(int k=0;k<n;k++){
					bii=bii+Y[i][k].a*U[k].b+Y[i][k].b*U[k].a;
			}
			return bii;
	}
	/**
	 * 求雅可比矩阵中的元素Hij
	 * @param 参数i是一个int型的数（行号）
	 * @param 参数j是一个int型的数（列号）
	 * @return 返回值是一个double型的数（雅可比矩阵中的元素Hij）
	 */
	public static double Hij(int i,int j){
			double Hij=0;
			if(i!=j){
					Hij=-Y[i][j].b*U[i].a+Y[i][j].a*U[i].b;
			}
			else{
					Hij=-Y[i][j].b*U[i].a+Y[i][j].a*U[i].b+bii(i);
			}

			return Hij;
	}
	/**
	 * 求雅可比矩阵中的元素Nij
	 * @param 参数i是一个int型的数（行号）
	 * @param 参数j是一个int型的数（列号）
	 * @return 返回值是一个double型的数（雅可比矩阵中的元素Nij）
	 */
	public static double Nij(int i,int j){
			double Nij=0;
			if(i!=j){
					Nij=Y[i][j].a*U[i].a+Y[i][j].b*U[i].b;
			}
			else{
					Nij=Y[i][j].a*U[i].a+Y[i][j].b*U[i].b+aii(i);
			}
			return Nij;
	}
	/**
	 * 求雅可比矩阵中的元素JRij
	 * @param 参数i是一个int型的数（行号）
	 * @param 参数j是一个int型的数（列号）
	 * @return 返回值是一个double型的数（雅可比矩阵中的元素Jij或Rij）
	 */
	public static double JRij(int i,int j){
			double JRij=0;
			if(i!=j){
					if(PQU[i].mark == 0){//PQ节点
					//if(i < s){
							JRij=-Y[i][j].a*U[i].a-Y[i][j].b*U[i].b;  //PQ节点时的计算公式
					}
					else if(PQU[i].mark == 2){//PV节点
					//else{
							JRij=0;  //PV节点时的计算公式
					}
			}
			else{
					if(PQU[i].mark == 0){
					//if(i < s){
							JRij=-Y[i][j].a*U[i].a-Y[i][j].b*U[i].b+aii(i);  //PQ节点时的计算公式
					} 
					else if(PQU[i].mark == 2){//PV节点
					//else{
							JRij=2*U[i].b;   //PV节点时的计算公式
					}
			}
			return JRij;
	}
	/**
	 * 求雅可比矩阵中的元素LSij
	 * @param 参数i是一个int型的数（行号）
	 * @param 参数j是一个int型的数（列号）
	 * @return 返回值是一个double型的数（雅可比矩阵中的元素Lij或Sij）
	 */
	public static double LSij(int i,int j){
			double LSij=0;
			if(i!=j){
					if(PQU[i].mark == 0){ //PQ节点
					//if(i < s){
							LSij=-Y[i][j].b*U[i].a+Y[i][j].a*U[i].b;  //PQ节点时的计算公式
					}
					else if(PQU[i].mark == 2){//PV节点
					//else{
							LSij=0;  //PV节点时的计算公式
					}
			}
			else{
					if(PQU[i].mark == 0){//PQ节点
					//if(i < s){
							LSij=-Y[i][j].b*U[i].a+Y[i][j].a*U[i].b-bii(i);  //PQ节点时的计算公式
					}
					else if(PQU[i].mark == 2){//PV节点
					//else{
							LSij=2*U[i].a;  //PV节点时的计算公式
					}
			}
			return LSij;
	}
	/**
	 * 求注入有功的不平衡量
	 * @param 参数i是一个int型的数（节点号）
	 * @return  返回值是一个double型的数（注入有功的不平衡量）
	 */
	public static double dP(int i){
			double dP=0;
			for (int j=0;j<n;j++){
			dP+=U[i].a*(Y[i][j].a*U[j].a-Y[i][j].b*U[j].b)+U[i].b*(Y[i][j].a*U[j].b+Y[i][j].b*U[j].a);
		}
		dP=PQU[i].p-dP;
		return dP;
	}
	/**
	 * 求注入无功或节点电压平方的不平衡量
	 * @param 参数i是一个int型的数（节点号）
	 * @return  返回值是一个double型的数（注入无功或节点电压平方的不平衡量）
	 */
	public static double dQU2(int i){
		double dQU2=0;
		if(PQU[i].mark == 0){//PQ节点
		//if(i < s){
			for (int j=0;j<n;j++){         //PQ节点时的计算公式
				dQU2+=U[i].b*(Y[i][j].a*U[j].a-Y[i][j].b*U[j].b)-U[i].a*(Y[i][j].a*U[j].b+Y[i][j].b*U[j].a);
			}
			dQU2=PQU[i].qv-dQU2;
		}
		else if(PQU[i].mark == 2){//PV节点
		//else {
			dQU2=PQU[i].qv*PQU[i].qv-U[i].a*U[i].a-U[i].b*U[i].b;  //PV节点时的计算公式
		}
		return dQU2;
	}
	/**
	 * 求雅可比矩阵
	 * @return 返回值是一个double型的二维数组（雅可比矩阵）
	 */
	public static double[][] J(){
		double[][] J=new double[2*n-2][2*n-2];
		for(int i=0;i<s;i++){     //求平衡节点号之前的各节点的2x2阶子阵
			for(int j=0;j<s;j++){
				if(Y[i][j].a!=0||Y[i][j].b!=0){
				    J[i*2][j*2]=Hij(i,j);
				    J[i*2][j*2+1]=Nij(i,j);
					J[i*2+1][j*2]=JRij(i,j);
					J[i*2+1][j*2+1]=LSij(i,j);
				 }
				else{
					J[i*2][j*2]=0;
					J[i*2][j*2+1]=0;
					J[i*2+1][j*2]=0;
					J[i*2+1][j*2+1]=0;
				}
			}
			for(int j=s+1;j<n-1;j++){
				if(Y[i][j].a!=0||Y[i][j].b!=0){
				    J[i*2][j*2-2]=Hij(i,j);
				    J[i*2][j*2-1]=Nij(i,j);
			        J[i*2+1][j*2-2]=JRij(i,j);
				    J[i*2+1][j*2-1]=LSij(i,j);
				}
				else{
					J[i*2][j*2-2]=0;
					J[i*2][j*2-1]=0;
					J[i*2+1][j*2-2]=0;
					J[i*2+1][j*2-1]=0;
				}
			}
		}
		for(int i=s+1;i<n;i++){     //求平衡节点号之后的各节点的2x2阶子阵
			for(int j=0;j<s;j++){
				if(Y[i][j].a!=0||Y[i][j].b!=0){
				    J[i*2-2][j*2]=Hij(i,j);
				    J[i*2-2][j*2+1]=Nij(i,j);
				    J[i*2-1][j*2]=JRij(i,j);
				    J[i*2-1][j*2+1]=LSij(i,j);
				}
				else{
					J[i*2-2][j*2]=0;
					J[i*2-2][j*2+1]=0;
					J[i*2-1][j*2]=0;
					J[i*2-1][j*2+1]=0;
				}
			}
			for(int j=s+1;j<n;j++){
				if(Y[i][j].a!=0||Y[i][j].b!=0){
				    J[i*2-2][j*2-2]=Hij(i,j);
				    J[i*2-2][j*2-1]=Nij(i,j);
				    J[i*2-1][j*2-2]=JRij(i,j);
					J[i*2-1][j*2-1]=LSij(i,j);
				}
				else{
					J[i*2-2][j*2-2]=0;
					J[i*2-2][j*2-1]=0;
					J[i*2-1][j*2-2]=0;
					J[i*2-1][j*2-1]=0;

				}
			}
		}
		return J;
	}
	/**
	 * 求不平衡量矩阵
	 * @return 返回值是一个double型的一维数组（不平衡量）
	 */
	public static double[] dPQU(){
		double[] dPQU=new double[2*n-2];
		for(int i=0;i<s;i++){
			dPQU[2*i]=dP(i);      //求注入有功的不平衡量   
			dPQU[2*i+1]=dQU2(i);  //求注入无功或节点电压平方的不平衡量
		}
		for(int i=s+1;i<n;i++){
			dPQU[2*i-2]=dP(i);
			dPQU[2*i-1]=dQU2(i);
		}
		return dPQU;
	}
	/**
	 * 求各节点电压修正量
	 * @return 返回值是一个double型的一维数组（各节点电压修正量）
	 */
	public static double[] dU(){
		double[] dU=Mat.multi(Mat.inv(J),dPQU);
		return dU;
	}
	/**
	 * 求各节点电压新值
	 * @return 返回值是一个Plural型的一维数组（各节点电压新值量）
	 */
	public static Plural[] New_U(){
		//double[] dU=dU();
		for(int i=0;i<s;i++){
			U[i].a+=dU[2*i+1];
			U[i].b+=dU[2*i];                  
		}
		for(int i=s+1;i<n;i++){
			U[i].a+=dU[2*i-1];
			U[i].b+=dU[2*i-2]; 
		}
		return U;
	}
	/**
	 * 判断收敛
	 * @param 参数aa是一个double型的数（收敛精度）
	 * @return 返回值是一个int型的数（即电压是否收敛，1表示收敛，0表示不收敛）
	 */
    public static int convergence(double aa){
    	int result=1;
    	//double[] dU=dU();
    	for(int i=0;i<2*n-2;i++){
    		if(Math.abs(dU[i])>aa){    //求各节点电压修正量的绝对值
    			result=0;              //若某个节点电压修正量的绝对值大于收敛精度aa
    			break;                 //则置result为0，既不收敛，强制退出循环
    		}	
    	}
    	return result;
    }
}

